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Abstract 

We study the relationship between chaotic behavior and the Central Limit Theo- 
rem (CLT) in the Kuramoto model. We calculate sums of angles at equidistant times 
along deterministic trajectories of single oscillators and we show that, when chaos is 
sufficiently strong , the Pdfs of the sums tend to a Gaussian, consistently with the 
standard CLT. On the other hand, when the system is at the "edge of chaos" (i.e. 
in a regime with vanishing Lyapunov exponents), robust (7-Gaussian-like attractors 
naturally emerge, consistently with recently proved generalizations of the CLT. 
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1 Introduction 



The relationship between weak chaos and the Central Limit Theorem (CLT) 
has been recently explored by several authors, both in low-dimensional chaotic 
maps [1,2,3,4] and in long-range many-body systems [5,6]. The Standard Cen- 
tral Limit Theorem (CLT) states that the (rescaled) sum of n independent 
random variables with finite variance has a Gaussian distribution in the limit 
n ^ oo. On the other hand, the theorem does not hold when the variables 
under consideration are weakly mixing and strongly correlated, as it happens 
in systems with zero or vanishing Lyapunov exponents, i.e. systems at the 
"edge of chaos". In order to consider this case, it has recently been intro- 
duced a q-generalization of the standard CLT, which states that for certain 
classes of strongly correlated random variables the (rescaled) sum approaches 
a g-Gaussian distribution [7]. The latter is a generalization of a Gaussian 
distribution which emerges in the context of the nonextensive statistical me- 
chanics [8]. It is defined as Gq{P,x) = A(q, — (1 — q)Px'^]^^^^^'^\ where q is 
the entropic index (such that for = 1 one recovers the usual Gaussian), /3 a 
parameter related to the width of the distribution and A{q, /3) is a normalizing 
constant. 

In previous works [5,6] we showed that g-Gaussian attractors exist for the 
rescaled sums of the angular velocities, calculated along deterministic trajec- 
tories of single oscillators (time-average Pdfs), in the weakly chaotic quasi- 
stationary (QSS) regime of the Hamiltonian Mean Field (HMF) model [9]. On 
the other hand, we also pointed out [10,11] several analogies ( met ast ability, 
continuous phase transition, chaotic behavior) between the HMF model and 
the Kuramoto model [12], probably due to their common origin from a more 
general model of coupled driven-damped pendula. In this paper we present 
new numerical results which extend these analogies to the CLT features. In 
particular we will show that also in the context of the Kuramoto model a 
very weakly chaotic microscopic dynamics with vanishing Lyapunov expo- 
nents gives rise to the emergence of robust g-Gaussian-like attractors, which 
disappear when the level of chaos increases, thus restoring the standard CLT 
behavior. 

In the first part of the paper we describe the Kuramoto model phase transition 
scenario, with specific attention to the chaotic behavior and to the role played 
by the distribution of the natural frequencies (uniform or Gaussian). In the 
second part we discuss the numerical results on the CLT behavior, which are 
strongly infiuenced by the interplay between chaos and synchronization. 
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2 Phase transitions, Chaos and MetastabiUty in the Kuramoto 
model 



The Kuramoto model [12,13,14] describes the following dissipative dynamics 
of a system of N sinusoidally coupled oscillators: 

K ^ 

Oi^oj, + -Y.sm{e^-e,), (1) 

i=i 



where 6i and ooi are, respectively, the angle (phase) and the natural frequency 
of the i-th oscillator, with i = 1, N, while i^T > is the coupling strenght. 
The natural frequencies are time-independent and randomly chosen from a 
symmetric, unimodal distribution g{u). They represent a driving term that, 
considering the oscillators as particles moving on the unit circle, forces each 
particle to rotate at its own frequency regardless of the other particles. On the 
other hand, the sinusoidal term, tuned by the coupling constant K, tends to 
synchronize the particles, so that, for values of K lower than a critical threshold 
Kc, the oscillators runs independently and incoherently, while for K > Kq 
the system undergoes a spontaneous phase transition towards a partially or 
fully synchronized state (depending on the value of K). Kuramoto [12] showed 
analitically that, for a given unimodal distribution g{uj), the critical value of 
the coupling follows the expression Kc = being g{0) the central value of 
the distribution. The degree of synchronization of the system is measured by 
the modulus < r < 1 of the dynamical vector order parameter r = re"^ = 
jf SjLi e*^^ , where (j) represents the centroid of the angles of the oscillators. 
Using the order parameter r, eq. (1) can be rewritten in its mean field version 
as 9i = LJi + Krsm{(f) — 9i), i = I, N, which allows to appreciate the role 
of both the natural frequencies uji and the absolute value of the product Kr 
on determining the synchronization degree of the Kuramoto stationary states 
(i.e. those with 9i = 0). In fact, as long as K is smaller than Kc, one has 
> Kr for all i and the oscillators drift at their own frequencies {drifting 
oscillators) along the unit circle, so that the system is incoherent and r ~ 0. 
When K increases above Kc, the system starts to synchronize and a cluster 
of locked oscillators, those with \uji\ < Kr, appears. In correspondence, r 
starts to grow (indicating partial synchronization) and the number of drifting 
oscillators decreases more and more until, for large value of K and for r > 0.8, 
the system becomes fully synchronized and all the oscillators rotate at the 
average frequency fl, a quantity which is conserved during the dynamical 
evolution. 

In ref.[ll], confirming and extending previous studies of other authors [15], 
we showed that the distribution g{u) plays a fundamental role for fixing the 
features of the phase transition in the Kuramoto model, even for large sys- 
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Fig. 1. (Lower panel) The asymptotic order parameter r of the Kuramoto model 
is plotted as a function of the coupling K for a system of = 20000 oscillators 
and for several bounded Gaussian distributions ^(w), with different standard de- 
viations a and u) G [—2,2]. Increasing a from 1 to 5 (triangles lines, from left to 
right), the 5'(w)'s changes continuously from a Gaussian-like to a uniform- like dis- 
tribution and, in correspondence, the phase transition changes continuously from 
second-order-like to first-order-like, while the critical value of the coupling Kc shift 
to the left. The curves corresponding to a true Gaussian distribution (open circles) 
and a true uniform one (open squares) are also reported for comparison. (Upper 
panel) The Largest Lyapunov Exponent is plotted as function of K for the same 
cases showed in the upper panel. Increasing a its behavior changes continuously 
(from left to right) from the characteristic, slowly decreasing, peak of a Gaussian 
g{ui) to the abruptly vanishing peak characteristic of the uniform one. Each dot 
represents an average over 10 runs. In the inset we show the power-law scaling of 
the LLE for small values of K. The LLE goes to zero as 7 • 10^ • K^'^. See text. 
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tern sizes. In particular we showed that the shape of the phase transition can 
be slowly changed from a continuous 2nd-ordcr-like one (very similar to the 
HMF phase transition and corresponding to a Gaussian g{u>)), to an abrupt 
Ist-order-like one (with coexistence of phases and corresponding to a uniform 
g{uj)), by tuning the standard deviation of a bounded Gaussian distribution 
of the natural frequencies. 

We report an example of such a behavior in Fig.l, where the phase transitions 
for a system oi N = 20000 are plotted for several values of the standard devia- 
tion a of bounded Gaussian distributions g{ujys, with to e [—2, 2] (examples of 
these distributions can be found in the insets of Fig. 2). In the lower panel, the 
asymptotic order parameter r is plotted as function of K for five increasing 
values of a (namely a = 1,1.5,2,3,5, see triangles lines). The limiting cases 
(7 = 1 and (7 = 5 approximate, respectively, the true Gaussian and the true 
uniform g{u) distributions, which are also reported for comparison (lines with 
open circles and open squares, respectively). 

In the upper panel, the behavior of the corresponding largest Lyapunov ex- 
ponent [LLE) is reported for all these distributions. As we already showed 
in ref.[ll], the partially synchronized regime of the Kuramoto model, char- 
acterized by the simultaneous presence of drifting and locked oscillators and 
by an order parameter < r < 0.8, exhibits a positive LLE with a peak 
just in coincidence of the phase transition, indicating the presence of strong 
chaos (see also [16]). On the other hand, in the incoherent phase, when all the 
oscillators arc drifting, we observe a vanishing (but positive) LLE, indicating 
a regime at the "edge of chaos". The inset of the upper panel of Fig.l reveals 
also an asymptotic power-law scahng of the LLE, which goes to zero as K'^ '^, 
a further analogy with a similar behavior observed in the HMF model [17], 
although in that case the power-law exponent was different. For the discussion 
of the next section it is important to note that, in the partially synchronized 
regime, the microscopic chaotic dynamics of the drifting oscillators enters in 
competition with the synchronization term represented by the locked oscilla- 
tors. The drifting term is prevalent provided that r < 0.2: we call this region 
of the partially synchronized regime the "fully chaotic" one, since it is also 
roughly characterized by a LLE > 0.01. On the contrary, for r > 0.2 the 
synchronization term becomes relevant and, even if the system is still sligthly 
chaotic (with a < LLE < 0.01), the cluster of locked oscillators starts to 
dominate the dynamics. Finally, for r > 0.8, we are in the fully synchronized 
phase and all the oscillators belong to the locked cluster which rotates at the 
average constant frequency Q (in both the panels of Fig.l, we report horizontal 
dotted lines in order to mark these different regions). Notice that in the fully 
synchronized regime the system is integrable and the LLE vanishes again as 
in the incoherent phase, but in this case it can also be slightly negative. Fur- 
thermore, being g{uj) symmetric with uj G [—2, 2], one has in this case ~ 0, 
therefore the locked cluster tends asymptotically to stop (we remind in this 
respect that the Kuramoto system is dissipative). 
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Fig. 2. Phase Diagram K versus a for the Kuramoto model with N = 20000. The 
critical line separates the fully synchronized phase (characterized by r > 0.8 and 
a slightly negative LLE) from the "edge of chaos" incoherent one (characterized 
by r = and a vanishing, but positive, LLE). The partially synchronized regime 
(characterized by < r < 0.8 and a positive LLE) is also visible between the two, 
around the phase transition. In particular, just below the critical line, we recognize 
a weakly synchronized sub-region of the partially synchronized regime, with values 
< r < 0.2 and LLE > 0.01, that we call "fully chaotic regime" (see text). In the 
insets, the bounded g{uj) distributions (with co € [—2,2]) used in the simulations 
are plotted for three increasing values of a. See text. 

A schematic picture of the previous scenario is plotted in Fig. 2, where the 
diagram K versus a is reported together with the labels of the various phases 
and the corresponding values of r and LLE (this picture extends the diagram 
presented in ref.[ll]). For a = 1 and a = 5 we have, respectively, the two 
estreme cases of Gaussian-like and uniform-like g{uj) (see the insets). The 
dashed critical line, drawn using the LLE peaks in Fig.l, corresponds to a 
2nd-order-like transition which - strictly speaking - would tend to a Ist-order- 
like one only in the limit a = oo, even if, as shown in Fig.l, the features 
of an abrupt Ist-order transition can be clearly recognized already for a = 
5. The partially synchronized regime is visible between the incoherent and 
the fully synchronized phases, together with its sub-region corresponding to 
the fully chaotic regime, which is recognizable immediately below the critical 
line. Finally, just above the critical line, a thick full line (blue in color) going 
across the partially synchronized regime indicates the presence of metastable 
states, see ref.[10]. In such states the interplay between locked and drifting 
oscillators hinders synchronization and does not allow r to reach immediately 
its maximum (stationary) value indicated by the curves in the lower panel of 
Fig.l. The situation is again similar to the one observed below the transition 
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in the quasi-stationary regime of the HMF model, where some particles are 
locked in a single cluster while the others move inchoerently on the unit circle 
(see ref.[ll]). 

In the next section we will expore the CLT behavior in the two limiting cases 
of a uniform and a Gaussian g{Lu) distribution, where true Ist-order and 2nd- 
order phase transitions occur. As shown in Fig.l, these two limiting situations 
are quite well approximated by the cr = 5 and a = 1 bounded Gaussians 
cases, but do not exactly coincide with them. In particular, for the uniform 
and Gaussian cases, the value of Kc is more consistent with the Kuramoto 
theoretical prediction, which becomes here, respectively, Kc = f ~ 2.54 and 

Kc — ~ 1.59, and both the r{K) and LLE(K) curves (open squares and 
open cirles lines) are slightly shifted on the right with respect to the a = 5 
and a = 1 ones. In the following we will refer to these curves to analyze the 
CLT features of the Kuramoto model in its three peculiar regimes: incoherent, 
fully chaotic and fully synchronized. 



3 CLT behavior in the Kuramoto model 

Stimulating by the analogies between the Kuramoto and the HMF models 
found in refs.[10,ll], and also by the observation of anomalous GLT behavior 

in the HMF context [5,6], we applied to the Kuramoto oscillators the procedure 
originally introduced in ref.fl] for constructing probability density functions 
(CLT Pdfs) of quantities j/'s expressed as a finite sum of stochastic variables, 
in order to identify possible violations of the standard CLT predictions. At 
variance with the HMF context, where we used angular velocities as stochastic 
variables, in this case the natural choice was to consider the angle variables 
of the N Kuramoto oscillators. Therefore we construct the rescaled sums y's 
by picking out, for each oscillator, n values of the angle 9i at fixed intervals of 
time 5 along the deterministic time evolution imposed by Eq.(l), i.e.: 

1 " 

Vi^^Y^UkS) with i^l,2,...,N. (2) 

k=i 

In such a way, the product S x n will give the total simulation time over 
which the sum is calculated. As explained in the introduction, for the stan- 
dard CLT, we would expect a Gaussian y's Pdf, but we already anticipated 
that the weakly chaotic microscopic dynamics, characterizing the Kuramoto 
model with vanishing LLE, will give rise to the emergence of g-Gaussian-like 
attractors, which disappear when the level of chaos increases. 
First we show the results obtained for a uniform g{uj) distribution of the nat- 
ural frequencies, then we will pass to consider a Gaussian one. In all the CLT 
numerical simulations we always consider a system oi N — 20000 oscillators. 
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Kuramoto Model- N=20000 - K=0.1 - Uniform g(a)) 
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Fig. 3. Uniform g{iv): CLT Pdfs for a Kuramoto system with = 20000 and K = Q.l 
(in this case Kc = 2.54), are reported for increasing values of 5 at fixed n (middle 
row) and for increasing values of n at fixed 5 (lower row): in both the cases robust 
fat tailed attractors appear, well fitted by g-Gaussians (full lines) with q = 1.7 and 
different values of /3. A standard Gaussian with unitary variance (dashed line) is 
reported for comparison in each plot. In the upper row of this and of the other 
figures by now on, the order parameter (on the left) and the LLE (on the right) are 
plotted for comparison as function of time, after a transient of 100 time-steps. See 
text. 

the same size used in Figs.l and 2, which offers a good compromise between 
statistical effectiveness and computational demand. We normalize also the y's 
to zero mean and unitary variance variables, for allowing a better comparison 
with standard Gaussian and g-Gaussian curves. 



3.1 Uniform g{u) distribution 



Let us start with the case of a uniform g{uj) distribution of the natural fre- 
quencies, with uj G [—2,2]. In the following we present several numerical sim- 
ulations, one for each region of interest in agreement with the phase diagram 
of Fig. 2 (see as an example the vertical section at a = 5, corresponding to the 
uniform-like g{uj) distribution case, quite similar to the true uniform case). 
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Kuramoto Model- N=20000 - K=2.53 - Uniform 
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Fig. 4. Uniform g{u;), N = 20000 and K = 2.53 (in this case Kc = 2.54). The fully 
chaotic regime tends produce Gaussian-like Pdfs, more visible increasing either 6 or 
n. Standard Gaussian curves are also reported for comparison. See text. 

In Fig. 3 we show the y's Pdfs of the same run a.t K = 0.1, i.e. in the "edge of 
chaos" incoherent phase. We plot the Pdfs for increasing values of 6 (200 and 
2000) at fixed ?7,=20 (middle row) and for increasing values of n (220 and 1000) 
at fixed 5=200 (lower row). The time evolution of both the order parameter r 
and the LLE is plotted in the upper row, after a transient of 100 time-steps. 
As expected, both of them stay around zero for the entire time evolution (in 
particular, LLE ~ 10~^). Correspondingly, one observes the occurrence of 
robust fat tailed Pdfs that persist increasing 6 or n. The latter can be fitted 
very well by g- Gaussian curves (full lines) with q = 1.7 and different values of 
(3 reported in the figure. A standard Gaussian with unitary variance (dashed 
line) is also plotted for comparison. 



Increasing K up to K = 2.53 (very near to the critical value Kc = 2.54 for the 
uniform g{uj) distribution and corresponding to the LLE peak in the upper 
panel of Fig.l - open squares line), one enters into the fully chaotic regime: 
here the system is still quite incoherent (0 < r < 0.2) but the high level 
of chaos {LLE > 0.01) drives the CLT Pdfs towards the standard Gaussian 
attractor, as shown in Fig.4. Again we plot both r{t) and LLE{t) in the top 
panels, while in the middle row and in the lower row we plot, respectively, 
the Pdfs of the same run for different values of 6 (200 and 3000) at n = 20 
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Kuramoto Model- N=20000 - K=3 - Uniform g(a)) 
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Fig. 5. Uniform g{io), N = 20000 and K = 3.0 (in this case Kc = 2.54). In the 
fully synchronized phase the CLT Pdfs coincide with the angles distributions of the 
oscillators in the locked cluster. Gaussian curves are also reported for comparison. 
See text. 

and n (220 and 1000) at 5 = 200. In this case, increasing either 5 or n, 
the Pdfs converge towards a Gaussian shape starting from the central part 
and slowly involving the tails. Again, a standard Gaussian is also plotted 
for comparison. It is important to stress that here the fully chaotic region fills 
almost all the partially synchronized regime and the Ist-order phase transition 
is characterized by a coexistence of the incoherent and synchronized phases 
which gives rise also (see full line in Fig. 2) to the appearance of metastable 
states [11] (in analogy with the QSS regime of the HMF model). In Fig. 4 we 
avoid these metastable states and we show a typical run whose LLE stationary 
value is consistent with the average value corresponding to the peak in Fig.l 
(i.e. LLE ~ 0.06 at K = 2.53). 



Finally, in Fig.5, we further increase the coupling up to X = 3, shifting in the 
fully synchronized regime. Here the locked cluster tends to include almost all 
the oscillators, whose angles do not change in time any more (we remind that 
the synchronized cluster tends to stop, being i7 ~ 0). In this situation one 
has r > 0.8 while the LLE is slightly negative, and the CLT Pdfs coincide 
(for large 5 or n) with rescaled angles distributions (being yi ~ 6i for each i = 
1, N). Summarizing, in the case of a uniform g{Lj) distribution of the natural 



10 



0,4 
0,2 


-0,2 
-0,4 





'"1 1 1 




oil 1 1 



T 



T 



T 



Kuramoto Model- N=20000 - K=0.6 - Gaussian g(a)) 
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Fig. 6. Gaussian g{u): CLT Pdfs for a Kuramoto system with = 20000 and 
K = 0.6 (in this case Kq = 1.59), are reported for increasing values of 5 at fixed 
n (middle row) and for increasing values of n at fixed 6 (lower row): as happened 
for the uniform case in Fig. 3, being in the "edge of chaos" regime, robust fat tailed 
attractors appear again, weh fitted by g-Gaussians (fuh lines) with q = 1.7. Standard 
Gaussians with unitary variance (dashed lines) are reported for comparison. In the 
upper row the order parameter (on the left) and the LLE (on the right) are plotted 
for comparison as function of time, again after a transient of 100 time-steps. See 
text. 

frequencies, the previous scenario confirms that only in the incoherent regime 
at the "edge of chaos", with vanishing LLE, robust g- Gaussian-like attractors 
appear in the CLT Pdfs of the Kuramoto model, then leaving room, when the 
level of chaos increases, to the standard Gaussian behavior predicted by the 
CLT. Let us to explore, now, what happens for a Gaussian g^u) distribution, 
when a 2nd-order phase transition occurs. 



3.2 Gaussian giu) distribution 



Considering a cr = 1 Gaussian g{uj) distribution of the natural frequencies, the 
region of interest in the diagram of Fig. 2 shifts towards the vertical section 
at 0" = 1 (corresponding to a Gaussian-like distribution, but quite similar to 
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Kuramoto Model- N=20000 - K=1 .65 - Gaussian g(co) 
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Fig. 7. Gaussian ^(w), N = 20000 and K = 1.65. The fully chaotic regime produces 
Gaussian-like Pdfs increasing either 5 or n. Gaussian curves are also reported for 
comparison. See text. 

the true Gaussian case). Notice that - at variance with what happened along 
the (T = 5 section - the partially synchronized regime is here mostly situated 
above the critical line, a footprint of the second order phase transition. As 
in the previous section we start studying the region of weak chaos in the 
incoherent phase. 



In the upper row of Fig. 6 we plot the time behavior of the order parameter 
and the LLE for a typical run a.t K = 0.6, and one can see that both vanish 
as expected (in particular LLE ~ 10~^). In the other panels CLT Pdfs are 
reported for increasing values of 6 (200 and 2000) at fixed ?t,=20 (middle row) 
and for increasing values of n (220 and 1000) at fixed 5=200 (lower row). As in 
Fig. 3, we found again robust fat tailed Pdfs that in this case can be fitted again 
very well by g-Gaussian curves (full lines) with q = 1.7 and slightly different 
values of (3. The usual standard Gaussian with unitary variance (dashed line) 
is reported for comparison. 

In the simulations showed in Fig. 7 we increase the coupling in the fully chaotic 
regime, up to the value K = 1.65, very near to critical point, where the mixing 
effects indicated by the peak in the LLE (see Fig.l, open circles line) slowly 
destroy the correlations and drive the CLT Pdfs towards a Gaussian-like shape 
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Kuramoto Model- N=20000 - K=4 - Gaussian g(co) 
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Fig. 8. Gaussian g{uj), N = 20000 and K = 4.0. In the fully synchronized phase the 
CLT Pdfs end up to coincide with the angles distributions of the oscillators in the 
locked synchronization cluster. Gaussian curves are also reported for comparison. 
See text. 



which again starts from the central part, as in the uniform g{co) case. Such 
a behavior is evident increasing either 6 {at n = 20) or n (at 5 = 200). In 
both cases we plot for comparison the usual Gaussian curve. In this situation 
no coexistence of different phases is observed, even if metastable states are 
still present for 1.65 < K < 1.75. Therefore, again, we avoid these states 
and choose a typical run whose LLE stationary value is consistent with the 
corresponding peak in Fig.l {LLE ~ 0.05 ai K = 1.65). 

Finally, in Fig. 8 we further increase K in the fully synchronized regime up to 
K = A and we observe a behavior analogous to that of Fig.5: being r > 0.8 
and LLE ~ but slightly negative, the CLT Pdfs tend to coincide with the 
angles distributions, assuming here a slightly asymmetric Gaussian-like shape 
for large values of 5 and n. Summarizing, also in the case of a Gaussian g{uj) 
distribution of the natural frequencies, only the vanishing Lyapunov exponent 
regime characterizing the incoherent phase gives rise to robust (/-Gaussian-like 
attractors in the CLT Pdfs. 
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4 Conclusions 



We have studied tlie relationsliip between CLT and cliaos in the Kuramoto 
model. By investigating, along the deterministic time evolution, the behavior 
of the sums of the phases of the oscillators, we have shown that when chaos is 
fully developed the CLT is satisfied and Gaussian Pdfs emerge as attractors. 
This is the expected result predicted by the standard CLT and obtained in 
other Hamiltonian models recently explored [5,6], but also in dissipative sys- 
tems like the logistic map [1,2]. On the other hand, when the system is at the 
"edge of chaos", i.e. well below the critical value of the coupling, where the 
LLE tends to zero, the standard CLT no longer applies. As observed in the 
recent literature for other models, in this case the Pdfs have very fat tails and 
the curves tend to attractors which can be very well fitted with q-Gaussians. 
Although we have presented only numerical simulations which cannot be con- 
sidered a rigorous proof, the results are in line with the recently generalized 
version of the CLT proposed by Tsallis [8] and provide further support to the 
claim that q- Gaussian-like attractors appear very often in situations where one 
has vanishing Lyapunov exponents so that mixing and ergodicity are violated. 
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